Center motions of nonoverlapping condensates coupled by long-range dipolar 
interaction in bilayer and multilayer stacks 
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We investigate the effect of anisotropic and long-range dipole-dipole interaction (DDI) on the 
center motions of nonoverlapping Bose-Einstein condensates (BEC) in bilayer and multilayer stacks. 
In the bilayer, it is shown analytically that while DDI plays no role in the in-phase modes of center 
motions of condensates, out-of-phase mode frequency (w ) depends crucially on the strength of 
DDI (ad). At the small-ad limit, u%(ad) — ^o(O) oc ad- In the multilayer stack, transverse modes 
associated with center motions of coupled condensates are found to be optical phonon like. At the 
long- wavelength limit, phonon velocity is proportional to y/a d . 

PACS numbers: 03.75.Hh, 03.65.-w 
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I. INTRODUCTION 

In contrast to the isotropic s-wave contact interaction, 
dipole-dipole interaction (DDI) has a very distinct char- 
acter - anisotropic and long-range. This character re- 
sults in various interesting phenomena in ultracold dipo- 
lar atom or molecule systems. The studies of DDI effect 
on the structure and dynamics of quantum many-body 
systems is at the forefront of both theoretical and ex- 
perimental interests. In recent years dipolar chromium 
( 52 Cr) atoms were created with a magnetic moment of 
6/ib {fJ-B is Bohr magneton), which is equivalent to a 
dipole moment d « 0.056Z? (ID ~ 3.335 x lCT 30 C-m) [J. 
More recently, creations of dipolar molecules have also 
been achieved in ultracold heteronuclear molecules such 
as 40 K 87 Rb with a much stronger dipole moment d « 
0.61? @, Q. From the theoretical point of view, many 
interesting physical properties have been studied, such 
as biconcave shapes of the ground state 041]) intriguing 
collapse mechanisms [5 , I7H1Q] , vortex structures [Tl| - [l3j , 
anisotropic solitons [14 1 and so on. 

Among many others, properties of nonoverlapping 
dipolar Bose-Einstein condensates (BEC) in bilayer or 
multilayer (quasi-lD optical lattice) stacks have at- 
tracted great attention in recent years. The topics under 
investigation include phonon instability fl5l. Hal , soliton- 
soliton scattering [lTfl 
ment condensation II 



pair superfluidity [181 ] . and hla- 
In these systems, condensates 
are effectively nonoverlapping between neighboring layers 
and as a matter of fact, contact or even short-range in- 
teraction plays no role between neighboring layers. How- 
ever, long-range DDI will still play an important role 
across different layers. Moreover in a real system, dipoles 
of the dipolar gas are usually aligned (by applying strong 
electric or magnetic field) along the direction of the 
stacks. This in turn will make dipolar molecular systems 
more stable against recombination and result in better 
opportunity of forming dipolar molecular BEC. 

This paper attempts to study collective excitations of 
nonoverlapping atomic or molecular BEC in bilayer and 
multilayer stacks, that are coupled by long-range DDI. 
The focus will be placed on the motions of the centers 



of condensates. When weakly interacting Bose conden- 
sates are loaded into a single harmonic (magnetic) trap, 
it is well known that, for small oscillations, center motion 
of condensates will undergo harmonic oscillation of fre- 
quency equal to the characteristic trap frequency [20l.l2lj]. 
This remains true even when the system has a long-range 
DDI in it. When condensates are loaded into a double 
well and form a bilayer system in, say, z direction (while 
a single harmonic trap is applied in the x-y plane) , center 
motions of the two condensate clouds can be very inter- 
esting upon the activation of long-range DDI. It will be 
shown explicitly that when DDI is present, both in-phase 
and out-of-phase motions of the two condensate centers 
will still undergo harmonic oscillations (for small oscilla- 
tions). However, while in-phase modes are independent 
of DDI, out-of-phase modes will depend crucially on the 
strength of DDI. 

In addition, it is also interesting to investigate center 
motions of nonoverlapping condensates in a multilayer 
stack (quasi-lD optical lattice along z-direction). By 
properly treating the boundary effect (see later), trans- 
verse modes associated with motions of the centers of 
condensates in each layer are found to be optical phonon 
like. At infinitely long wavelength limit (q z — ¥ 0), phonon 
mode frequency just equals to the harmonic trap fre- 
quency in the x-y plane and at the long-wavelength limit 
(q z zo -C 1, #o is the lattice constant), phonon mode ve- 
locity is found to be proportional to the square root of 
the strength of DDI. Due to the long-range character of 
DDI, for the multilayer system of finite number of lay- 
ers, it is important to treat properly the boundary effect. 
With this regard, we have introduced a truncation num- 
ber (Nj!) which corresponds to number of neighboring 
layers included for satisfactorily converged results. It is 
found that Njl will depend on the ratio of lattice constant 
(zq) and condensate radius in the x-y plane only. 

To end this introductory section, we emphasize two 
things. Firstly, the results presented in this paper, es- 
pecially the dependence of DDI strength on center mo- 
tion mode frequencies, are believed to be valid even when 
the system is not Bose condensed. Since long-range DDI 
survives in the non-condensed system, it is encouraging 
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that the experiment can also be done on the ultracold 
heteronuclear dipolar molecules which exhibits a strong 
dipole moment but nevertheless is yet to be Bose con- 
densed @, 111 ■ Secondly, the strength of DDI can actu- 
ally be extracted through measurements of out-of-phase 
mode in the bilayer or transverse phonon mode in the 
multilayer stack. 

The paper is organized as the following. In Sec. HH 
energy functional of a nonoverlapping bilayer system is 
given. Proper trial wave functions are introduced within 
the variational framework and minimized ground-state 
energies are obtained. In Sec. IIII1 analytical results of 
the mode frequency for center motions of condensates are 
given. It is shown that while in-phase mode frequency is 
independent of the strength of DDI, out-of-phase modes 
depend crucially on the strength of DDI. In Sec. IIV1 we 
extend the study to the center motions of condensates in 
a multilayer stack. By properly treating the boundary 
effect, transverse phonon modes are obtained. At the 
long- wavelength limit, we show that phonon velocity is 
proportional to the square root of the strength of DDI. 
Sec. IV] is a conclusion. 



II. NONOVERLAPPING BILAYER SYSTEM 

We consider a nonoverlapping bilayer system with 
same BEC atoms or molecules in each layer. Energy 
functional of the system is given by 



E — Ei + Ei + Ei 



where (i — 1, 2) 



Ei = N J drVf(r) 
N - 1 



2m 



V 2 + T4xt,i(r) 
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and 



ff |^(r)| 2 + / dr'V dd (v - r')|^(r')f 



(1) 



(2) 



E 12 =N 2 / / drdr'V dd (r-r')\Mr)\ 2 \Mr')\ ■ (3) 



Here Ei and E 2 correspond to intralayer energies and 
E12 is the interlayer energy caused by DDI. N is the 
number of condensed atoms or molecules in each layer 
and tpi is the normalized wave function for layer i. 
g = 4:TTh 2 a/m with a the s-wave scattering length and 
V dd (r) = d 2 (l — 3 cos 2 #)/|r| 3 is the dipole-dipole interac- 
tion with 8 the angle between r and the dipole orienta- 
tion. For magnetic dipoles, the strength d 2 = fion^/Air 
with fi m the magnetic dipole moment, while for electric 
dipoles, d 2 — d 2 /47re with d e the electric dipole mo- 
ment. For a nonoverlapping bilayer considered in current 
context, condensate wave functions are not overlapped 
across the two layers, i.e., J dr|V'i(r)| 2 |?/'2(r)| 2 ~ 0. As 
a consequence, interlayer energies associated with s-wave 
interaction as well as interlayer hopping are neglected. 



To create a nonoverlapping bilayer ultracold BEC sys- 
tem, one can set up a double well with a large bar- 
rier in the middle. One possible candidate of nonover- 
lapping double potential well in z direction is V — 
\m {ijJ 2 x 2 + uj 2 y 2 + (jJ 2 z z 2 ) + Vi cos 2 (irz/di) by applying 
a harmonic trap together with a deep optical (lattice) 
trap [22| | . Alternatively, one can apply a harmonic trap 
together with a repulsive barrier potential to obtain 
V = \m(Lo 2 x x 2 +uj 2 y y 2 +Lu 2 z z 2 ) + V exp(-z 2 /£ 2 z ) 
By controlling the values of Vi vs. de or V vs. t z , 
effective nonoverlapping bilayer system can be formed. 
As a matter of fact, wave function of each layer is ex- 
tremely narrow (pancake like) in z-direction and to the 
leading order, trap potentials experienced by the atoms 
or molecules in each layer can be approximated by the 
following anisotropic harmonic potentials, 



Vc 



ext,l,2 



(r) 



f 



LU 2 x x 2 +LU 2 y 2 +LU 2 (z T zi/2) 2 ■ (4) 



Here zi, determined by the experimental setup, corre- 
sponds to the spacing between the two trap minima lo- 
cated at z = ±zi/2. For nonoverlapping bilayer systems, 
lo z is typically much larger than oj x and uj v . In this paper, 
as depicted in Fig. [TJ we consider the dipolar system to 
which dipoles are oriented along z direction. With this 
geometry, the system will be the most stable compared 
to others. 

Within the variational framework, Gaussian ansatz is 
used for the trial wave functions of the bilayer system. 
As long as the system is not close to a collapsed state, 
Gaussian function should be a good trial wave function 
for the system [24|. For layer 1 and 2, we then take 



^1,2 (r) = Aexp 




,(5) 



where variational parameters R x , R y , and R z corre- 
spond to condensate radii in x, y, and z directions, 
A = 1/ \J R x R y R z ir 3 / 2 is the normalization constant, and 
Zq corresponds to the distance between the two conden- 
sate centers located at z — ±zq/2. Note in general that zq 
can be different from zi [the latter corresponds to trap 
minima, see Eq. Q]. When DDI vanishes, zq will be 
identical to z\. However, when DDI is present, zq will be 
smaller (larger) than zi if DDI is attractive (positive) in 
z direction. Nevertheless, for the present nonoverlapping 
bilayer system, zq should be not much different from z\. 

Truly speaking, due to anisotropic nature of DDI, real 
wave function for each layer will not be symmetric in z 
direction even though a symmetric potential trap ((3]) is 
in effect. However, when uj z 3> LJ x ,uj y , wave function in 
each layer is still relatively symmetric in z direction. Sub- 
stituting trial wave functions §5§ and trap potentials (0| 
into energy functional ([I])-©, we obtain 
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FIG. 1. Sketch of the nonoverlapping bilayer condensates 
and their center motions along the a;-direction. Red arrows 
denote the direction of dipoles. Blue stacks represent the 
stationary condensates, while black dash lines represent their 
movements. Frame (a) corresponds to the in-phase mode and 
frame (b) corresponds to the out-of-phase mode. 
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In ([6]) and throughout this paper, we have rescaled the 
energy E/(huj x ) — > E and the time tui x — > t. Be- 
sides, all the lengthes (R x , R y , R z , z , z\) are scaled by 
the magnetic length in cc-direction, I = yjhj muj x , and 
the ratios X y = u y /uj x and A 2 = lj z /lj x . Dimension- 
less coupling strength a d = d 2 m/h 2 £ and a s = aji. 
F = F(zo, R x , R y , R z ) is related to the inter layer energy 
due to DDI, while F = F(z = 0,R x ,R y ,R z ) is related 
to the intralayer energy due to DDI. More explicitly, for 
dipoles aligned along the z direction, F is given by the 
following integral 



F = -^ I </k, x P 



1 



2 v x 



cos ( k z — 



3k 2 



klRl/Rl + k 2 y R 2 /Rl + fcj 



(7) 



Values of R x , R y , R z , and zq for the ground state are 
obtained by solving the conditions dE/dRi = (i = 
x 7 y, z) and dE/dz — 0. 



III. CENTER MOTIONS OF BILAYER SYSTEM 



By variational approach, suitable dynamical variables 
should be added into the trial wave functions (|5|). In 
Ref . [13, HH , it was shown in the single-trap system that 
equations of motion for the center of condensates are de- 
coupled from equations of motion for the width of con- 
densates. For a bilayer system studied in this context, 
similar decoupling occurs although the proof is omitted 
for brevity. As a matter of fact, for the motions of the 
center of condensates, dynamical wave functions can be 
taken to be 



ipl t2 (r,t) = Aexp 



V 2 {z T V2) 2 



2R 2 



x exp ■ 



[x - 2:1,2ft)] 
2R 2 . 



2R 2 

2 



ixd(t)}, (8) 



where dynamical variables are added associated with the 
center motions only. Here Xi(t) corresponds to the fluctu- 
ation of the center of condensate of layer i in x-direction, 
while Ci (t) corresponds to the sloping phases of the con- 
densates of layer i. The task is to find the equations of 
motion and solve it. 

We start from the effective Lagrangian, L = T + E, 
where T is given by 



J .7=1,2 V 7 



dt 



dt 



(9) 



and E is given by Eqs. (d])-©- Substituting Eq. © into 
the Lagrangian and expanding dynamical variables up to 
second order (for small oscillations), one obtains 



L = 



E 

i=l,2 



N 



EL 
2 



N 2 a d G{ Xl -x 2 ) 
2R x R y R z 



where G = G(zq, R Xl R y , R z ) and given by the integral 



G 



1 

6^2 



dkki 



exp 



1 - 



2 v 



k 2 



3k 2 



cos 



k z ^ 

R z 



\R?JRI + k 2 y R 2 /R 2 y + k 2 



(11) 



It is interesting to note that s-wave scattering coupling 
a s has completely dropped out in L in (|10|) . Equations 
of motion can be derived using the Lagrange equation, 
It %i = §§ ' wnere 1 can b e an y one °f t ne f° ur dynam- 
ical variables. By assuming q = qoexp(iu)t), after some 
algebra we obtain two branches of excitation modes: 



As mentioned before, in this paper we focus on the mo- 
tions of the center of condensates along x direction (see 
Fig. [TJ. Generalization to considering as well the center 
motions of condensates along z direction is straightfor- 
ward. However, for condensates strongly confined in z 
direction, condensate wave functions are rigid in z di- 
rection and as a consequence collective excitation in z 
direction will cost more energy. Besides observation of 
these motions of much smaller amplitude are relatively 
more difficult. 



<4 - 1, 

lo 2 = 1 + Na d 



2G 



(12) 



where Wj corresponds to the in-phase mode associated 
with the solutions of X\ = x 2 and c\ = c 2 [such as the 
motion depicted in Fig. [TJa)] . While oj corresponds to 
the out-of-phase mode associated with the solutions of 
x\ = —x 2 and ci = — c 2 [such as the motion depicted in 
Fig. [lb)]. 
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FIG. 2. The square of out-phase mode, w„, plotted as the 
function of Nad (N is particle number in each layer and 
ad is DDI coupling) for A z = 16, 49, and 100 respectively. 
The insets show linear dependence of Nad at small Nad- 
In frame (a), s-wave coupling JVa s = and in frame (b), 
Na a — 10. 2i is chosen to be 1.21 in both frames. 



Since x-direction harmonic trap frequency uj x is the 
energy unit used in this context, the results in (|12l) show 
that in-phase mode frequency of x-direction center mo- 
tions of bilayer condensates is just lo x , independent of the 
long-range DDI. For the out-of-phase mode, in contrast, 
the mode frequency is shifted from uj x and the deviation 
depends crucially on the strength of DDI (ad)- To be 
more explicitly, the deviation will depend on a d , the num- 
ber N, as well as the four lengths (R x ,Ry,R z ,zo). Since 
number N and the four lengths are measurable quanti- 
ties, one can then calculate G [using Eq. (fTTj)] with the 
measured values of N and (R x ,R y ,R z ,Zo). This simply 
means that once center motions of the bilayer is per- 
formed and out-of-phase mode frequency lu is measured, 
one can actually extract the value of ad via Eq. (fl~2"j). 

Fig. [2] shows the square of out-of-phase mode ui 2 as 
the function of Nad- For convenience, we consider the 
isotropic case in xy plane (R x — R y = R and \ y = 1), 
and have chosen z\ — 1.21. Fig. Uta) plots the Na s = 
case with X z chosen to be 16, 49, and 100 respectively. 
While Fig.[2][b) plots Na s = 10 case with the same choice 
of the three A z 's. The insets show the linear dependence 
of Nad when Nad is small. Within variational frame- 
work, values of R, R z , and zq are determined by min- 
imizing the energy functional ([5]). It is useful to check 
that when both s-wave and DDI couplings are absent, 
a s = ad = 0, one finds that zq = z\, R = 11, and 
R z = (1/4)4 (1/7)4 an d (1/ 10 K corresponding respec- 
tively to X z — 16, 49 and 100 cases. However, when ad 
is present, zq becomes shorter than z\, while R and R z 
become larger than those for the case of ad = 0. This sim- 
ply means that the ratio of zq / R z will become smaller in 
the presence of ad- Nevertheless, in the strong confin- 
ing regime (z n /R z > 4) considered in this context, the 
change of zq/ R z due to ad is minor. For example, in case 
of a s = and Nad = 30, ratio z n /R z is found to be 4.07, 
7.9, and 11.6 for X z = 16,49 and 100, as compared to 
4.8, 8.4, and 12 for the case of a s = ad = 0. 

When results of Fig. &b) are compared to those in 



Fig.[UJa), one sees that repulsive s-wave coupling a s acts 
to reduce the deviation of ui from lu x . While short-range 
a s plays no role between neighboring layers, its repulsion 
actually increases the radii of condensate (R x , R y , R z ) in 
each layer and consequently Gj R x R y R z becomes smaller. 
This indicates that ui 2 or its slope against Nad (at small 
Nad) will be smaller. In addition, it is found that regard- 
less of the value of a s , for the same Nad, the larger X z is, 
the closer tu is to w x . In 52 Cr atom dipolar BEC, it has 
been measured that d 2 m/h 2 ~ 24 A. If atom number in 
one layer is N ~ 10 4 and the harmonic oscillator length 
I ~ 1 /im, then it is estimated that Nad ~ 20. This gives 
a reference how large u> is when the results of Fig. [2] are 
considered. 



IV. NONOVERLAPPING MULTILYAER STACK 

In this section, we extend to study center motions of 
condensates in a multilayer stack. As before, condensates 
are assumed to be non-overlapping between neighboring 
layers. The kind of system can be realized in a deep 
quasi-lD optical lattice 25}. In an analogous way, energy 
functional of the multilayer system can be given by 



E 



E K+E £ » 



where 

E n = N I ,lr, ■■(!•) 

N -1 



-^V 2 + I4 xt (r) 
Am 



9\ipn(r)\ 



dr'V dd (r-r')\i> n (r' 



A|2 



(13) 

(14) 
ipn(r) 



N 2 / / drdr'^(r-r')|^„(r)| 2 |^ n+m (r')|1;i5) 



and 

E„ , 



Here N s corresponds to the number of layers in the stack 
and N corresponds to the number of atoms in each layer. 
E n represents the intralayer energy for layer n, while 
E n m represents the interlayer energy between layers n 
and (n+m) coupled by DDI. The external trap, V cx t, con- 
sists of magnetic and optical traps, where magnetic trap 
is m(uj 2 x 2 +oj 2 y 2 ) /2, while optical trap is sE r sin 2 (7rz/zo) 
with zq the spacing between neighboring layers, s the 
strength of optical trap, and E r = Ti 2 n 2 /2mz 2 the recoil 
energy. 

In the calculation, we shall use the following approxi- 
mation for E in Eq. ([TUl) : 



E 



E„ 



(16) 



|m| — 1 



where N c is a truncation number to which how many 
neighboring sites are included for E n<m (\m\ = 1 corre- 
spond to the two nearest neighbors). With (IT5|) . periodic 
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FIG. 3. Sketch of the multilayer stack and the motions of the 
center of condensates. Frame (a) corresponds to a q z — in- 
phase phonon mode and frame (b) corresponds toai} 2 = tv/zq 
out-of-phase phonon mode. Red arrows denote the direction 
of dipoles. Blue stacks represent the stationary condensates 
and black dash lines represent the motions of condensates. 
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FIG. 4. Square of the transverse phonon modes plotted as 
the function of (kzo/n) 2 . Two cases, (Na B ,Nad) = (0,1) 
and (0.5, 1) are presented. The inset shows the linearity at 
long wavelength limit (kzo <C 7r) and the slope is equal to the 
square of phonon velocity. 



boundary condition (PBC) which connects the two ends 
is applied. Owing to the long-range nature of DDI, it is 
expected that for a satisfactory converged result, N c ^> 1. 
On the other hand, for a trusty result, it is required that 
the truncation number should be much smaller than the 
total number of layers in the stack (N c <C N s ). 

In the current multilayer stack, for simplicity, we also 
apply the Gaussian ansatz for the trial wave function 
of each layer [similar to that of Eq. ((5J]. It should be 
emphasized again that due to the boundary effect, true 
wave function associated with each layer is not perfectly 
symmetric in z direction. Nevertheless, for the current 
nonoverlapping condensates under study, a symmetric 
wave function will be a good approximation for each 
layer. After a lengthy derivation, we obtain the following 
energy functional for the multilayer system 



E_ _ f J_ J_ t J_ sE, 



+ — + — + — (1 - e ~^ R ^ /x h 

~ A D2 X A D2 X O K 1 C J 



R 



H 2 

x , x 2 "'y 



N - 1 



4 " 4 2R x R y R z 




a s — adF 



N 2 a d 
2R x R y R 



where 



■r(iv c ), 



r(^ c ) 



E 



F„ 



(17) 



(18) 



I' 



with F m being just F in Eq. {7} with cos(k z zo/R z ) term 
replaced by cos(mk z zo/ R z ). Once again, it is assumed 
that dipoles of the dipolar gas are aligned along the z- 
direction and we consider only the motions of the center 
of condensates in x direction for small oscillations (see 
Fig. [3]). In ([T7|> . all energies and lengths are rescaled in 
the same way as those in the bilayer case. The recoil 



energy will then reduce to E r = it 2 £ 2 /2zq in the dimen- 
sionlcss unit. 

With the same variational approach, we obtain the La- 
grangian of the system (keeping dynamical variables up 
to second order) 



l =ny: 



■t'n 

2 



ENa d (x 
. „o „ „ Lrr7 



'■1 = 1 



4R x R y R z 



(19) 



with G m being just G in Eq. (fTTj) with cos(k z zo/R z ) 
term replaced by cos(mk z z /R z ). Similar to those in 
Eq. ([5]), dynamical variable x n = x n (t) corresponds to 
the fluctuation of the center of condensate of layer n in 
x-direction, while c„ = c n (t) corresponds to the slop- 
ing phases of the condensates of layer n. Assuming that 
In = qoexp(iu)t — nkzo) (q n represents any one of the 
dynamical variables in layer n and k z —> k for brevity 
afterwards), we obtain the following dispersion relations 
for the transverse modes: 



w 2 (fc) 



1 



Na r , 



R x RyR 



■A(iV c ), 



where 



A(AQ 



E 

|m|=l 



[1 — cos(mfczo)] Gr. 



(20) 



(21) 



The above mode is analogous to an optical phonon mode 
in crystals to which cj(k — 0) corresponds to a in-phase 
mode, while ui(k = tt/zq) corresponds to an out-phase 
mode for neighboring layers. These two special cases are 
illustrated in Fig. EI It can be simply checked that u>(k = 
0) = 1. That is, in-phase mode frequency is just lj x . 
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Moreover, if we include only the nearest-neighbors (N c = 
1) for A in Eq. (|2ip . out-of-phase mode uj(k = tt/z q ) will 
be exactly the same as lo for the bilayer system [see 
Eq. (|12[)] . At the long- wavelength limit (kzo 1), one 
obtains 



uj 2 (k) ~l + v 2 (N c )k 2 
where the square of phonon velocity 

Na d z 2 



v 2 (N c ) 



t^j> -f^y ^~^z 



E 

|m| — 1 



777 G Tl 



(22) 



(23) 



Thus v 2 oc Nad at the low-fc limit. Measurements of 
dispersion relations of the phonon modes thus can give 
direct information on the value of DDI. 

In Fig. [U we plot the dispersion relations of transverse 
phonon mode uj with (Na s , Nad) = (0,1) and (0.5,1) 
respectively. The numbers used are based on assuming 
that number of 52 Cr atoms in each layer is N — 400, 
magnetic length £ is about 1/xm, and hence Nad is about 
1. The curves presented in Fig. H] are obtained using 
a proper truncation number N c leading to satisfactory 
converged results (see later). Moreover, we assume that 
sE r = 300 and zo = 0.7£ and hence s is about 30, which 
is in the deep optical lattice regime. By minimizing the 
energy functional (fTT|) . we obtain that R z /R ~ 0.074 
and zq/R z ~ 7.0 for the (Na s , Nad) = (0,1) case and 
R z /R ~ 0.061 and z /R z ~ 7.0 for the (Na s ,Na d ) = 
(0.5, 1) case. Inset of Fig. @] shows the linear dependence 
of k 2 on uj 2 at the long wavelength limit. The slope is 
equal to the square of phonon velocity, v 2 (N c ). When 
the curve of finite a s is compared to that of a s — 0, one 
sees that repulsive s-wave coupling a s acts to suppress the 
phonon mode as well as the phonon velocity. While short- 
range a s plays no role between neighboring layers, its 
repulsion actually increases the radii of condensate (R x , 
Ry, R z ) in each layer and consequently K(N C )/ R^.R y R z 
becomes smaller. This indicates that w 2 as well as the 
phonon velocity will be smaller. 

Finally the behaviors of T(N C ), A(N C ), and v 2 (N c ) as 
the function of N c are studied. It is important to first 
note that r(AT c ), A(N C ), and v 2 (N c ) all exhibit the same 
converging behavior. In Fig. EJa), T(N C ) is plotted as 
the function of N c for fixed R z — 0.1£ and R = 1£ and 
three choices of z = 0.5£, 0.9£, and 1.3& It is seen 
clearly in Fig. [SJa) that T(N C ) converges at some value 
of N c and the larger the zq is, the smaller the N c is for 
the converging result. To be more explicit, we define a 
critical value AT? of N r such that 



T(Ni)-T(Ni-l) 



< 10 



-3 



(24) 



T(N?) 

In fact, Af? = 19 and 25 respectively for the two curves 
in Fig. |H 

Fig. Hth) plots AT? as a function of z /R. It is found 
that N§ depends only on the ratio of zq/R regardless 
of the value of R z . This occurs because for the current 
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FIG. 5. (a) r(N c ) plotted as the function of N c for three 
choices of zo = 0.5£, 0.9£, and 1.3£ and fixed R — 1£ and 
R z — 0.1£. (b) Critical number N§ (see text for definition) 
plotted as the function of zq/R. 



nonoverlapping multilayer system, R z -C zo and R z is no 
longer a well-defined length scale owing to the long-range 
character of DDI. As a matter of fact, A^ depends only 
on the ratio of Zo / R- When the ratio z /R is larger, the 
corresponding N§ is smaller and A^ — s- 1 in the limit of 
large zq/R. The results of N§ indeed can help clarify- 
ing whether our results of transverse phonon modes are 
trusty or not. How does it work? Let us consider in a 
real experiment that total number of layers is N s = 100 
and the ratio Zq/R is about 1. According to the results of 
Fig-EIb), N-? = 16 for z /R = 1. In this case, one does 
meets the criterion N c <§; N s and our results of trans- 
verse phonon modes are trusty when a real experimental 
measurement is compared to. 



CONCLUSIONS 



In this paper, analytical solutions of mode frequencies 
for the motions of the center of condensates are studied in 
bilayer and multilayer (quasi-lD optical lattice) stacks. 
In the bilayer, it is shown that while DDI plays no role in 
the in-phase modes of center motions of condensates, out- 
of-phase modes (u) ) depend crucially on the strength of 
DDI (ad)' More explicitly, lo will depend on condensate 
radii of each layer (R x , R y , R z ), interlayer spacing (zq), 
as well as Nad- Therefore one can actually extract the 
value of a d if w Q , (R Xl R y , R z ), and zq are measured ex- 
perimentally. In the multilayer stack system, transverse 
(optical) phonon modes and phonon velocity are derived 
explicitly. Proper treatment was made for the bound- 
ary effect and it turns out that the truncation number 
to which how many neighboring sites should be included 
for the long-range DDI is a function of zq/R (R is the 
condensate radius in the transverse direction) only. 
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